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ABSTRACT: 

The  equations  governing  free  turbulent  mixing  are  derived  from  the 
Navier -Stokes  equations  and  transformed  into  a  mathematical  plane  which 
is  explicitly  independent  of  the  eddy  viscosity  model.  The  coupled 
momentum  and  turbulent  kinetic  energy  equations  are  analytically  solved 
in  the  transformed  plane  by  a  perturbation  technique  and  subsequently 
retransformed  into  physical  space  based  on  a  hypothesized  dependence  of 
the  eddy  viscosity  on  the  turbulent  kinetic  energy.  The  adequacy  of  a 
given  model  in  reproducing  the  velocity  and  turbulent  kinetic  energy 
field  is  assessed  by  comparing  the  results  of  the  analysis  with  some 
experimental  data  of  planar  turbulent  wake  mixing  in  constant  adverse 
and  favorable  pressure  gradients. 
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NOMENCLATURE 

A  =  Cross-Sectional  Area  of  Wind  Tunnel  Test  Section 

p  =  Strength  of  Pressure  Gradient,  Equations    (U9)   and   (51) 

C  =  Constant   in  Equation   (32) 

Dv  =  Turbulent  Dissipation 

6*  =  Displacement  Thickness 

o 

e  =  Perturbation  Parameter  =  l/u   ;  DK/p  in  Equation  2k 

F  =  Arbitrary  Function  in  Equation  (73) 

F*  =  Fe-Jf^)dcP 

G  =  Greens  Function 

H  =  Shape  Factor  =  §*/e 

I  =  Initial  Conditions   in  Equation   (75) 

i  =  Unit  Vector  in  x-direction 

J  =  Nonhomogeneous  Function  in  Equation   (7^)   =  g(cp>Y)e  J    ^'^P 

j  =  Unit  Vector  in  y-direction         


/      2  2  2\ 


K  =  Turbulent  Kinetic  Energy  =  \ 

k  =  Constant  in  Equation   (52) 

k  =  Unit  Vector  in  z-direction 

2       -2 

X  =  u       -  u 

*  e 

I  =  Macros cale  of  Turbulence 

(j,  =  Coefficient  of  Viscosity 

[i,  =  Eddy  Viscosity 

v  =  Kinematic  Viscosity 

v,  =  Turbulent  Kinematic  Viscosity 

P  =  Pressure 

Y  =  Stream  Function  Defined  in  Equation   (25) 

cp  =  Transformed  Independent  Variable   in  Equation   (33) 

q  =  Dynamic  Pressure 


IV 


p  =  Density- 
ax,  =  Turbulent  Diffusion  Coefficient  in  Equation  (23) 
K 

o 

S  =  Velocity  Ratio  Function  =  (u  /u  ) 

t  =  Time 

t  =  Apparent  Turbulent  Shear  Stress 

9  =  Momentum  Thickness 

u  =  Axial  Velocity  (x-direction) 

U  =  Average  Axial  Velocity  External  to  the  Wake 

v  =  Lateral  Velocity  (y-direction) 

V  =  Velocity  Vector  =  ui  +  vj  +  wk 

w  =  Transverse  Velocity  (z-direction) 

x  =  Axial  Coordinate 

y  =  Lateral  Coordinate 

z  =  Transverse  Coordinate 

Subscripts 

e  =  External  to  the  Wake 

i,«j,k  =  Three  Orthogonal  Components  of  a  Vector 

L  =  Exit  of  Test  Section 

0  =  Inlet  of  Test  Section 

t  =  Turbulent  Function 

0,1  =  Zeroth  and  First  Order  Solutions  of  x  and  K 

Superscripts 

u  =  Time  Average  of  u 

u'  =  Fluctuating  Part  of  u 


INTRODUCTION 

The  solution  of  the  Navier-Stokes  equations  as  applied  to  the  problem 
of  turbulent  flow  has  classically  been  approached  either  from  the  point  of 
view  of  elegant  rigor  on  limitingly  simple  flows  or  one  of  empirically 
guided  analysis  of  flows  of  real  interest.  The  Reynolds  averaging 
technique  is  nearly  a  pre -requisite  to  attacking  any  real  turbulent  flow 
and  the  associated  necessity  to  empirically  close  the  system  of  equations 
with  a  Reynolds  stress  model  makes  the  mathematics  tractable  yet  often 
unreliable.  Various  so-called  eddy  viscosities  which  arise  from  the 
Boussinesq  laminarization  of  the  apparent  turbulent  stresses  are  functions 
of  the  flow  field  and  as  yet  there  are  no  known  universal  functions  which 
adequately  model  these  stresses  in  all  cases. 

Possibly  the  most  successful  application  of  eddy  viscosity  approaches 
lies  in  the  area  of  free  mixing  where  velocity  differences  through  the  flow 
field  are  small  and  the  associated  turbulent  field  is  a  fairly  simple  one. 
However,  in  the  pressure  gradient  situations  which  are  encountered  in 
ejector  and  combustor  mixing  phenomena,  classical  eddy  viscosities  fail 
since  they  are  explicitly  independent  of  the  turbulent  field  which,  in 
these  cases,  can  become  complex  and  exert  a  dominant  effect  on  the  form  of 
the  eddy  viscosity.   Modern  approaches  have  attempted  to  include  the  local 
turbulence  structure  effect  on  the  eddy  viscosity  through  an  explicit 
dependence  on  the  local  turbulent  kinetic  energy. 

In  conjunction  with  a  highly  idealized  wake  mixing  experiment,  the 
two-dimensional  incompressible  turbulent  wakes  in  constant  adverse  and 
favorable  pressure  gradients  have  been  studied  analytically  with  a  formu- 
lation of  the  equations  which  allows  the  coordinated  solutions  of  the 
velocity  and  turbulence  field  in  a  transformed  plane  which  is  explicitly 


independent  of  the  eddy  viscosity  model.   Subsequently,  the  retransform- 
ation  of  the  equations  for  a  variety  of  models  allows  for  a  comparison 
of  the  adequacy  of  the  models  and  an  evaluation  as  to  which  most 
accurately  reproduces  both  the  velocity  and  kinetic  energy  fields. 

EQUATION  DEVELOPMENT 
The  basis  for  any  rigorous  analysis  of  problems  involving  turbulent 
fluid  flow  must,  to  the  best  of  current  knowledge,  be  founded  in  the 
Navier -Stokes  equations.   For  an  incompressible,  Newtonian  fluid  these 
may  be  written 

7  •  V  =  0  (1) 


P  ^  =  -  VP  +  n  V^  (2) 


No  analytical  solution  of  the  full  equations  appears  possible  and, 
although  some  success  has  been  shown  with  numerical  approaches  to  the 
solution  of  the  equations,  as  of  yet  no  numerical  attack  for  a  genuinely 
turbulent  flow  is  feasible  (Ref.  l) .  Apart  from  the  numerical  approach, 
only  Fourier  analysis  has  shown  any  progress  in  the  solution  of  the 
equations.  However,  with  this  method,  only  limitingly  simple  turbulent 
fields  have  been  treated  and  its  relevance  to  a  general  turbulent  mixing 
problem  has  yet  to  be  proven. 

Classically,  the  most  fruitful  approach  in  the  analysis  of  turbulence 
has  been  to  decompose  each  of  the  dependent  variables  in  Equations  1  and 
2  into  a  mean  term,  which  is  independent  of  time  or  has  a  long  character- 
istic period  with  respect  to  the  turbulent  fluctuations,  plus  a  fluctua- 
ting term  whose  time  average  is  zero.   The  velocity  field 

V  =  ui  +  v j  +  wk  (3) 


may  be  decomposed  term  by  term  such  as  the  decomposition  of  u 

u  (x,y,z,t)  =  u  (x,y,z)  +  u'  (x,y,z,t)  (k) 

When  this   decomposition   is   applied  to  each  term  in  the   continuity 
equation  we  obtain 


dji  +  £l  +  ^  +  QL  +  *lL  +  *lL  =  o  (5) 

dx      dy      dz        dx        dy         dz 


If  we  now  take  the  time  average  of  Equation   (5),  we  obtain 


|H+|v+|w=0  (6) 

dx       By       dz 


which  is  the  continuity  equation  for  the  mean  flow.   Upon  subtracting 
Equation  (6)  from  Equation  (5),  we  obtain 


*i£  +  dw  +  ^-°  ™ 


This   is  the  continuity  equation  which  the  velocity  fluctuations  must 
satisfy.      Prior  to  applying  this  decomposition  technique  to  the  momentum 
equations  we  first  reformulate  the   convective  operator  in  a  conservative 
form,  with  the  aid  of  the  continuity  equation.     The  general  convective 
operator 


i  +  ul  +  vl  +  wl  (8) 

St  dx  dy  dz 


becomes 


&♦&»<>+&*)♦£«(   )  O) 


where  the  appropriate  dependent  variable  is  placed  within  the  parentheses, 
Applying  this  formulation  to  Equations  (2)  and  expanding  them  in  rectang- 
ular Cartesian  coordinates,  we  obtain 


3u        3     /   2x         b     /      s        9     /      >  1  dp  2 

tt  +  tt  (u   )   +  r-  (uv)   +  —  (uw)   =   -  -  ^  +  v  v  u 
3t        dx  By  3z  p    3x 


g^w^('?)^w--|t^  (10> 


dw^d/N         d     /      n         B/2x        1    Sp  2 

dT  +  a^  (uw)  +  a7  (vw)  +  a^  (w  }  =  p"  dT  +  v  v  w 

Each  of  the  dependent  variables   in  Equations    (10)    is  now  decomposed  and 
the  appropriate  expression  is   inserted  into  the  equations.      The  resulting 
system  of  equations   is 


du  +  du'    +  3u     +       3uu'    +  3u'      +  3uv  +  3uv'    +  3u'v  +  3u'v' 
at         at  3x  3x  dx  3y  3y  3y  3y 


duw       3u*w       duw'        du'w'    _        1  dp        1  dp'    .  2-  2 

+    — +    — +    — +    — -    -    —  — - +vVU    +    vVU 

dz  dz  dz  dz  p   dx        p     dx 


—  —  —         —  —2—2 

d_v       d_v_|_  +  duv       du'v  +  duv'    +  du'v'    +  3v_  +  „  3w'    +  dv' 

at      at       ax       ax        ax         ax        ay  ay        ay 

(ii) 

d^  +  dwv  +  d^vi  +  aw^i  =  _  i  dF  _  i  a?:         2-         2  , 
az       az        az         az  p  ay     p  ay  v 


dw  ,  3w'   duw  J  du'w  ,  duw'  ,  du'w'  ,  dvw  JL  dvw'   dv'w  ,  dv'w' 
+  — — —  +  +  +  +  — _    +  +  +  + 


dt    dt    dx    dx     dx     dx     dy    dy     dy     dy 


dw2  dww'        dw'2  1  dp       1  dp'  „2-   .  2    , 

+    rr +    2    -r +    — = r-£ T"-  +vVW+vVw' 

dz  dz  dz  p   dz        p     dz         v  v 


If  we  now  take  the  time  average  of  each  equation,  we  can  simplify  the 
system  by  dropping  out  all  those  terms  whose  time  average  is  zero  (i.e. 
those  which  have  terms  that  are  linear  in  the  fluctuating  properties). 
Note,  however,  that  we  must  retain  the  nonlinear  terms  containing 
fluctuations  since  the  products  or  powers  of  purely  fluctuating  terms  may 
generate  a  steady  time  averaged  value.  When  the  appropriate  time  averages 
are  taken  of  each  term  in  Equations  (ll)  we  retain  the  following  system  of 
equations 


du    d  /— 2\    3  / — \    3  / — \     Id/—,    ,2\    d  — ; — r   3  — ; — r  ,    2— 

at  +  dx"  (u  }  +  d7  (uv)  aT  (uw)  =  "  p~  ax"  (p  +  pu  }  -  d7  u  v  "  dT  u  w  +  v  v  u 


dv    d  / — n    d  /-2x    d  / — x     Id/-,    ,2\         d  -7— r   d  T7  j   „2_ 

dt  +  ^  (uv)  +  a7  (v  }  +  di  (w)  =  "  p  d7  (p  pv  }  "  ^  u  v  "  dT  w  v  +  -  v  v 


If  +  I-  (w)  +  #-  (w)  +  I-  (w2)  =  -  -  |-  (p  +  pw'2)  -  |-  u^T7"  -  |-  v^w1"  +  v  V2w 


(12) 

Note  that  the  price  of  this  "simplification"  through  the  Reynolds  averaging 
has  been  the  introduction  of  six  new  unknowns  by  discarding  of  all  of  the 
phase  information  of  the  fluctuations.  For  a  steady,  two-dimensional  mean 
flow  Equations  (12)  may  be  simplified  to 


|H  +  ^  =  0  (13) 

dx   dy  v  J 


-  Bu    ,         Bu  13/—.         ,2v  B     — : — r    ,  2--  /_,  x 

u  to  +  v  b7  =  "  *  &  (p     pu    }  _37  u  v      v  7  u  (lU) 


-  Bv   ,   —  Bv  1   3    /r  .        i2n        B    — s — r  .         2-  /,_* 

u  —  +  v  —  = —  (p  +  pv'    )-—  u'v'   +  v  v  v  (15) 

Bx  By  p  By  Bx  v  v   " 

In  general  we  may  make  the  assumption  that  the  apparent  pressures    (pu'    , 

2  —  — 

pv'    )   may  be  neglected  with  respect  to  p.     We   should  note  that  when  p  = 

constant  some   care  must  be  exercised  in  applying  this   approximation  since 


-  —  — -  (pu'    )   and  -  —  —  (pv'    )   may  be   significant  terms   in  the  equations 
p  Bx  p  By 

for  a  particular  problem.      In  addition,   for  a  "thin"  wake,  we  may  make  the 
standard  boundary  layer-type  parallel  flow  approximation  that 

*-  A  dU 

|P  =  f£  =   _  pU         e  (   6) 

dx       $x  *   e    $x  v      ' 

This  eliminates  the  necessity  of  solving  the  lateral  momentum  equation, 
since  the  only  unknowns  are  u  and  v  with  p  being  imposed  on  the  mixing 
region  by  the  external  flow.   The  only  remaining  term  which  explicitly 


involves  the  turbulent  field  is  -u'v'  which  represents  an  apparent  Reynolds 
shearing  stress  (t). 


t  =  -  pu'v'  (17) 

In  general  this  apparent  stress  is  very  large  with  respect  to  the  average 
laminar  shear  stress  and  we  may  neglect  the  laminar  component.  With  these 
approximations  the  resulting  system  of  equations  to  be  solved  is 


|^  +  ^=0  (18) 

Bx   By 


Bx    By    p  dx   p  By 


6 


At  this  point,  some  empiricism  is  necessary  since  t  is  not  retrievable 
from  the  equation  set  we  have  derived.  This  is  clear  from  Equations  (ll) 
when  either  is  multiplied  by  the  appropriate  fluctuating  velocity  and 
time  averaging  is  performed  to  obtain  an  equation  for  the  Reynolds  stress, 
triple  correlations  appear  which  are  unknowns  also.  This  cascade  of 
unknown  higher  order  correlations  continues  for  all  further  equations 
which  are  derived.  Some  success  (Ref .  2)  has  been  made  attacking  these 
problems  by  making  purely  heuristic  approximations  in  higher  order 
correlations  in  an  effort  to  make  the  least  sensitive  approximations 
possible  in  the  equations.   However,  in  lieu  of  attacking  this  spiralling 
set  of  equations,  it  has  generally  proved  more  effective  to  introduce  some 
empirically  based  models  of  the  Reynolds  stress  into  Equations  (l8)  and 
(19).  The  most  successful  models  have  been  based  on  extensions  of 
Prandtl's  mixing  length  analysis  which  analogizes  turbulent  eddy  momentum 
transfer  with  the  molecular  manifestation  of  viscosity.   In  line  with  this 
approach  an  eddy  viscosity,  ^ , ,  is  introduced  into  the  problem  via  the 
definition 

H  '-  -  »^f  (20) 

where  (j,,  is  a  function  of  the  local  mean  flow  field.   Classical  models 
infer  that  u.,  is  purely  a  function  of  the  local  non-turbulent  mean  flow, 
however,  more  modern  results  (Ref.  3)  indicate  that  some  specific 
dependence  on  the  turbulence  field  is  indicated.  Thus  we  hypothesize  that 

M-+  =  M-+  (u>v,K) 


where 


K  =  I  (u'2  +  v'2  +  w'2)  (21) 


In  order  to  implement  a  model  of  this  form,  we  must  formulate  and  solve  an 
equation  for  K  along  with  the  momentum  and  continuity  equations.   In  order 
to  obtain  an  equation  for  K,  we  multiply  Equations  (ll)  by  u' ,  v1 ,  w' 
respectively,  add  them  and  take  the  time  average.   This  operation  results 
in  the  following  equation  for  K  which  is  most  conveniently  written  in 
tensor  notation 


Dt        ax.  Vp     3        2     k   i   3  i  Vax.      ax. ;; 


(22) 


,du.    du.s      au\   du*.v2 

2 1  (_i  +   <A  _   f i.  +  J 

dx .    9x./   v  \  ax .    dx. 


The  respective  convection,  diffusion,  production,  and  dissipation  terms 
have  been  modelled  by  Patankar  and  Spalding  (Ref .  4)  in  the  following 
equation 


r-V 


-  dK  +  -  dK  =  d_ 

dx    dy  "  ay  Lav   Sy_ 


t  dK" 

K 


-.2   D 


+  vt(|)  "f  <*> 


Where  av  is  the   effective  Prandtl  number  for  the  diffusion  of  K  and  D     is 
K.  ft 

the  turbulent  dissipation.     With  appropriate  auxiliary  expressions   or 

equations   for  u . ,  D    ,   and  a     we  may  consider  the  system  of  equations  to 
t         ft  ft 

be   closed.      Clearly  the  adequacy  of  the  equation  system  in  modelling  any 

particular  flow  field  is   dependent  on  the  exact  formulation  which   is  used 

for  the  unknown   coefficients . 

Based  on  dimensional  reasoning,  we   can  specify  the   coefficients  ^,. 

and  J)     to  be 
ft 

(j,,  =  const,  x  density  x  velocity  x  length 

D  =  const,  x  density  x  velocity  /length 

8 


From  experimental  results  (Ref.  3)  we  specify  a     to  be  a  purely  empirical 

constant  in  the  range  0.5  -•  1.0. 

With  the  explicit  expression  which  we  will  test  for  p,  ,  along  with 

a  value  for  a    ,   only  the  formulation  for  D„  is  left  unspecified.   Jones 
ft  Jv 

(Ref.  5)  has  hypothesized  an  equation  for  e  =  D  /p  °f  the  following  form 
t,     >   ,v+  a   v  du.  .du.   Qu     C0e 

Dt  '*  ax.  w  ax  J   i  k  vt  ax.  Vax.  ax./   k       (2k) 


3        e    j  j    1     j 


where  C-,  and  r     are  empirical  constants.   In  lieu  of  this  complication, 

we  will  make  some  experimentally  justified  approximations  for  D  and  its 

K 

dependence  on  K  and  the  mean  velocity  field. 


MATHEMATICAL  ANALYSIS 
I.  Equation  Reformulation 

In  order  to  reduce  Equations  (18) ,  (19)3  and  (23)  to  a  form  more 
amenable  to  an  analytic  approach,  we  introduce  the  stream  function 


Y,  =  u   ,  Y   =  -  v  (25) 


which  automatically  satisfies  the  continuity  equation.  With  this  definition 
of  the  stream  function,  the  von  Mises  transformation  is  applied  to  the 
independent  variables  (x,y)  to  transform  them  into  (x,y)  via  the  following 
equations 


k'h-*h  w 


(26) 


-^-  =  u  — 

ay  "  ay 


When  this  transformation  is  applied  to  Equations  (19)  and  (23),  we  obtain 


3uz.  _J_dp  +  _3_/   _3u 


dx     pu  dx   dY  \  t 


*  %)  ™ 


dK   3  /vtU  3K 


-.2   D 


+ 


Vj 


u  (P)     -  -i  (28) 


dx   dy  \  aR  BY/    t  VdY/    pu 

where  continuity  need  not  be  solved  since  v  does  not  appear  explicitly  in 
the  equations  when  u  and  K  are  expressed  in  stream  function  variables. 

The  dependent  variable  u  is  now  transformed  into  X  via  the  equation 

X  h  ue2  -  u2  (29) 


Upon  substituting  u  =  (u   -  x)2  into  Equations  (27)  and  (28)  we  can  write 
the  equations  of  momentum  and  turbulent  kinetic  energy  as 

&->*».  C1  -  ^  ft  (30) 

U  SY 


9x 


1  -Ion  -1 


% &DH1  - ^)'  If] +  ^M1  - ^)2 ft)  - -hC1  - 5)  1 <31) 


'K  u 


ku  U  LL.U  N  U 

e  e  e  ^t   e  e 


To  this  point,  based  on  the  assumptions  previously  outlined,  these 

equations  are  as  exact  as  the  specifications  of  Bvt   \>,,   and  a    .      To  be 

DK 
consistent  with  our  previous  dimensional  reasoning,  the  coefficient  —  in 

H 

the  kinetic  energy  equation  may  be  expressed  dimensionally  as 


DK  2 

—  =  constant  x  (velocity/ length) 


10 


Following  experimental^  indications  we  hypothesize  the  exact  formulation 
for  this  coefficient  to  be 


i  -   C  x  K//  (32) 


where  i   is  the  maxroscale  of  the  turbulent  flow,  generally  accepted  to  be 
of  the  order  of  the  mixing  width  and  a  physical  measure  of  those  turbulent 
eddies  most  intimately  involved  in  momentum  transfer.  We  have  not  yet 
explicitly  used  an  independent  expression  for  v  in  the  equations  and,  as 
we  shall  see,  the  hypothesis  for  its  formulation  need  not  be  specified 
until  the  equations  have  been  solved  in  a  transformed  plane. 

With  these  equations  and  expressions,  we  again  transform  the  indepen- 
dent variables,  this  time  from  (x,y)  to  (cp,Y)  using  the  following  defini- 
tion of  cp 

x 
cp  =   f  u  v.  dx  (33) 

d0  e  z 

where  we  have  introduced  a  crucial  yet  common  and  experimentally  justified 

assumption  that  y+  =  v+(x)  alone.  However,  our  interest  lies  in  specifying 

a  useful  v+(K)  which  is  implicitly  v+  [K(cp>y)]  in  "the  transformed  plane  and 

we  are  free  to  select  a  particular  y.  for  optimum  results,  thus  v+(x)  is  in 

J  t 

actuality  v^  =  vt  [K(ep,Y.)]. 

Applying  this  transformation  to  the  equations ,  we  can  write  the 
momentum  and  turbulent  kinetic  energy  equations  as 


ue    9y 
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dK     i   a 


i  _1  p  1 

2    am  1  /  v     \"2/Av\^  f!  /  v     \  "2 


H*tL\  ^)     3Yl       Uu2l  U2J     UJ  ,2U2^  2;  (35) 

e  e  e  e  e 


Subject  to   the  following  boundary  and  initial  conditions 
B.C.x->0;K->0     as     y  -  =  °° 


(36) 


I.C.   x   =  Xi(Y)      ;     K  =  K±(y)   at  cp  -  0 


II.  Equation  Simplification—Linearized  Case 

For  many  wake  or  jet -like  parallel  flows  the  velocity  within  the  wake 
differs  only  slightly  from  that  in  the  local  external  flow.   In  these  cases 


-*2  «  1 

u 
e 


and  we   can  write  Equations    (3^)   and   (35)   as 


&  =  d_JL 


^       a/ 


(37) 


2T 


dK  .      1    d^K       /_C_ \  v    ,       1      /3X\  /Qox 

^  K  ^  £   u  4u 


e  e 


subject  to  the   identical  boundary  and  initial  conditions   specified  for 
Equations    (3*0   and   (35). 

III.      Equation  Simplif ication—Nonlinearized  Case 

In  free  mixing  cases,   although  the   square  of  the  velocity  defect   (^) 

2 
is   small  with  respect  to  u      ,    it   is  not  negligible   and  some  influence   of 

the   finiteness   of  the  velocity  defect  must  be   included.     With  the  governing 
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equations  in  the  (cp,Y)  plane  it  is  particularly  convenient  to  include 
this  effect.  This  is  accomplished,  first  with  the  momentum  equation,  by- 
expanding  the  dependent  variables  in  powers  of  the  parameter  which  makes 
X  small,  namely  some  measure  of  the  freestream  velocity. 
With  the  following  definitions 


2 

B-(J)  ...Jg  (39) 

e        u^ 
0 


the  momentum  equation  becomes 


*  ■  (i  -  .  s  x)4  s-%  (to) 

y  dY 


We  are  seeking  here  only  a  first  order  correction  to  the  linearized  set  of 
equations,  however,  successive  higher  order  approximations  may  be  derived 
in  exactly  the  same  manner.   Using  the  binomial  theorem,  we  expand  the 
diffusive  term  and,  keeping  only  a  first  order  correction  to  the  linearized 
case,  we  can  write  the  momentum  equation 


£*  (i  - *.. -a  x)  j4  (*i) 

T  dY 

00 

We  now  expand  the  dependent  variable  x  i-n  a-  power  series  x  =  2  e  x  • 

i=0     X 
and,  in  line  with  the  objective  of  obtaining  a  first  order  correction  to 

the  linearized  system,  we  retain  only  the  first  two  terms  in  the  series 

which  results  in 


0   ®   1 
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Upon  substituting  this  expression  into  Equation  (hi)   and  collecting  terms 
of  equal  order  in  e  we  obtain  the  following  equations  for  yn»  Xn 


(te) 


ax0 

*2*o 

3cp 

.    2 
dY 

dcp 

.    2 
dY 

~   2   S  x0         p 

By 

(^3) 


Clearly  the  zeroth  order  equation  is  the  strictly  linearized  case  and  the 
first  order  solution  is  the  correction  term  for  the  finiteness  of  the 
velocity  defect  which  tends  to  force  the  solution  to  satisfy  the  full 
equation. 

Similarly,  the  linearized  solution  for  K  may  be  corrected  for  the 
improved  velocity  field  solution  and  for  the  direct  effect  of  a  finite 
velocity  defect.   Here  again  we  seek  a  first  order  correction  for  the 
linearized  solution  by  expanding  ^  and  K  in  Equation  (35)  in  terms  of  e . 
The  expressions  for  u  are  expanded  with  the  binomial  theorem  with  use  of 
the  appropriate  expressions  for  x  =  Xn  +  e  Xtj  the  appropriate  zeroth  and 
first  order  solutions  to  the  momentum  equation.   Clearly  the  equation  for 
K  is  now  explicitly  dependent  upon  e  and  we  can  expand  K  also  in  powers  of 
e,  again  to  first  order,  so  that  a  correction  to  the  linearized  case  may  be 
obtained.   Upon  substituting 


K  =  KQ  +  «  Kl 


the  respective  equations  for  K  and  K  can  be  obtained  by  collecting  terms 
of  equal  order  in  e .  The  equations  for  K  and  K  can  be  written 


1U 


2  2 

e        e 


2 
dK           ..     3   K.. 

*V                Oy-        ^2 

C        K              1 

A2  1    2ctk 

e 

,       1      dX0  SX1 

c     s  *o  Ko 

2u  2    dy      dY 
e 

,2     2           2 
e 

s  ~ 


a  r    o 


x 


By  Lao  5y  J 


2 

+  ^(itHif)   ('5) 


IV.  Specification  of  u  (cp) ,  i(cp) 

To  this  point,  whether  a  linearized  or  nonlinearized  approach  is 
taken,  the  solution  can  be  developed  analytically  only  providing  u  ,  I 
are  known  functions  of  cp.  We  want  to  generate  solutions  for  x(cp>Y)>  K(cp,y), 
hypothesize  a  dependence  of  v+  upon  K  and  unwind  the  transformation  via  the 
definition  of  cp  in  the  following  manner. 

dcp  =  v+  u     dx  (U6) 


dcp 


r  — nrt tt  =  f  tt    ax  (W 

J    v+    (K  [cp,Y,])        J      e  v      ' 


which  results  implicitly  in 

f(cp)  =  g(x) 

After  performing  this  integration,  we  can  solve  for  cp(x)  based  on  a  given 
u  (x)  and  the  hypothesized  functional  dependence  for  v  (K) .  The  solutions 
can  then  be  retransformed  into  physical  variables  and  the  validity  of  the 
hypotheses  can  be  checked  by  comparison  with  the  experimental  u  and  K  field, 
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In  order  to  obtain  expressions  for  u  and  I   as  functions  of  cp,  we 
must  first  specify  a  general  flow  field  to  be  studied.  The  following 
development  serves  only  as  a  guide  toward  the  specification  of  a  class  of 
functions  which  u  (cp)  and  jfc(cp)  must  be  in  order  to  adequately  model 
realistic  wake  flows  in  pressure  gradients.   However,  none  of  the  approxi- 
mations in  the  following  section  involves  an  actual  specification  of  v+ 
since  u  and  Z   could  be  specified  as  any  arbitrary  class  of  functions  of 
cp  and  the  actual  physical  flow  field  could  later  be  inferred. 

The  flow  external  to  the  wake  in  the  experimental  phase  of  the 
investigation  was  flowing  through  a  channel  whose  area  ratio  is 


f"  =     1  (W) 


0  J   1  +  px 


where 

2 


p  =  -^ —  (w 


For   incompressible  flow,   this   area  distribution  results   in  a  velocity 
distribution 


u 
e 


—  =  J  1  +  f3x  (50) 


U0 


with  a  resultant  pressure  gradient 


dp 


dx 


pqo  (51) 


In  order  to  map  u  and  I   into  the  cp  plane ,  we  need  a  reliable  approximate 
guide  as  to  the  shape  of  the  v+  function  so  as  to  specify  the  transformation 
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from  x  to  cp.   From  Schlichting  (Ref  6)  planar  incompressible  constant 
pressure  wakes  may  "be  approximately  solved  with  a  v+  °f  "the  fori 


•m 


vt  =  k  ue  0  (52) 


In  this  expression  8  serves  as  the  characteristic  length  and  for  weak 
wakes  is  a  measure  of  the  macroscale  (i)  as  is  shown  by  the  following 
development.  From  the  definition  of  y  we  can  write 


l 
e  1 


I-  -,       Y„  ,       x-2 


•-J'^.-J'C1--^  d*  (53) 


0   v    u 
e 


which  for  weak  wakes  can  he  approximated  by 


u 


^Je(i  +  -^x)dY  ») 


~0       2u 
e 


In  addition  from  the  definition  of  momentum  thickness 


Y 
0  ~e      "e/       0  "e      "e' 


•  ■irf^-D^-jI'f^-f)*       (") 


the  following  approximate  expression  can  be  derived 


Ye   1 


e     0  2u 


By  comparing  Equations  5^  and  56 


u  £  =  Y     +  u  9  (57) 

e  e         e 


Since  mass  flux  in  the  wake  and  hence  f  will  be  approximately  conserved, 
9  will  reflect  the  functional  dependence  of  the  macroscale  on  x. 
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In  order  to  obtain  an  expression  for  9  we  can  return  to  the  momentum 
integral  equation  written  with  the  specified  external  field  which 
generates  a  constant  pressure  gradient  —  =  -  Pqn.  The  momentum  integral 
equation 


dU 

i +  f  if  tH  +  s  =  °  <58) 

e 


can  be  integrated  for  a  constant  shape  factor  (H)  resulting  in 


60  =  (1  +  px)1  +  H/2  (59) 


With  this  external  velocity  field,  the  Schlichting  model  for  v  ,  and  the 
solution  for  9  we  can  write  the  transformation  as 

X   2 
cp  =  k  P  u   9  dx  (60) 

J0  e 

2 
after  inserting  the  expressions  for  9  and  u   we  obtain 


2 

ke^u 


cp  =  p 


°o  0   r/,   Q  a  -  h/2   .-]  ,,.  v 

(1  -H/2)L(1+PX)        -1]  (6l) 


Using  this  relationship,  the  expressions  for  u  (x)  and  9(x)  can  be 
transformed  into  u  (cp)  and  0  (cp) 

/  6  \X/(2    "  H) 

".W""^1  •  7    2   f      u2  0  (62) 


,  o  N(H  +  2)/(H   -   2) 

e(cp)  =  e0  (1  -       2  p 2  cp)  (63) 


[   -  2/  k9oU0 
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It  should  now  be  clear  that  we  can  consider  the  system  of  governing 
equations  in  the  problem  to  be  closed  and  well  defined  since  all 
functional  coefficients  have  been  specified  as  functions  of  cp.  Namely 
S(cp)  and  l(<p) 


n"  -  ®  ■ (1  -  <ste  •> 


2/(H   -  2) 


(Sh) 


and  from  Equation   (57) 


4(cp)ue(cp)    =  Ye  +  u090  (l   -  g    P  g  cp) 

\k^~2J  k0ouo 


(H  +  1)/(H   -  2) 


(65) 


The  specification  of  the  value  of  the   shape  factor    (H)  to  be  used 
must  be   consistent  with  the   formulation  and  the   inherent  approximations 
From  the  definition  of  H  we  can  write 


6* 

H  =  —  = 


J(i-f)*      rKi-t)* 


tHi-B*  J't(1-6 


dy 


(66) 


H  = 


J(l-§* 


(67) 


iC( 


1  - 


H  = 


_2L- 

u 
e 


-  1 


dY 


j[r-(i-^']-a, 


(68) 
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Expanding  the   integrand,   consistent  with  the  order  of  the  analysis  we 
obtain 


S[1  +  ^2-i 


2u 


dy 


H 


11-1  +  -*-=■  I  df 
L  2u     J 


=  1 


(69) 


If  we    insert  this  value   of  H  into  Equations    (6l),    (6U),    and   (65)   we   can 
write 


cp 


2kVo 


(1  +  px)2    -  1 


(70) 


B(«p)-(J)     =(l  + 


2kVo 


2  V 


(71) 


4(cp)  ue(cp)  =  Ye  +  u0e0  (1  + 


2ke0u0 


2  V 


-2 


(72) 


SOLUTION  OF  THE  EQUATIONS 
The  general  form  for  each  of  the  linearized  and  nonlinearized  equa- 
tions  for  both  K,  x   can  ^e  written 


U  =  ^|  +   f^)F   +   g(q>,   f ) 


(73) 


where  f(cp)  =  0  for  the  momentum  equations. 
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If  we  let  F*  =  Fe~J  ^P'^,  simple  differentiation  verifies  that  the 
solution  to  the  equation 

9   ay 

satisfies  the  general  form  of  Equation  (73). 

The  general  solution  of  Equation  (7*+)  can  be  written  in  terms  of 
Greens  function  G(x,y;cp,Y)  as 


+oo  CD       +00 

F*(cp,Y)  -  J      G(0,y;cp,Y)   I  (y)  dy  +  J     J      G(x,y;cp,Y)  J  (x,y)  dx  dy       (75) 

-00  0         -co 

a       a2 

where  the  Greens  function  appropriate  to  the  —  -  -^-—  operator  is 

p 

!»(»*».»)  -        p         «p  [-  l&^f]  (76) 

7  2rr(cp   -  x)  VVK 

Once  the  appropriate  initial  conditions  (i)  and  nonhomogeneous  functions  (J) 
are  inserted  and  integrated,  the  solution  of  F  is 


Jf(cp)dcp  (7?) 


F  =  F*e 

COMPARISON  WITH  EXPERIMENT 
In  order  to  test  the  adequacy  with  which  the  solutions  presented  in 
Equations  (73)  >  (75)  5  and  (77)  predict  a  turbulent  free  mixing  flow  field, 
a  sample  quadrature  was  performed  which  could  be  compared  with  experiment. 
Figure  1  presents  the  initial  conditions  for  x  and  K  which  were  used  in 
Equation  (75)  to  solve  for  xn>  X-,  >   K  ,  and  K  in  constant  adverse  and 
favorable  pressure  gradients.  Figures  1,  2,  and  k   contain  the  results  for 
a  constant  favorable  pressure  gradient  with  the  experimental  conditions  in 
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indicated  on  the  figures.  Figures  1,  3>  and  k  contain  the  results  for 
the  same  calculations  with  a  constant  adverse  pressure  gradient  of  the 
same  magnitude.  The  mathematical  retransformation  plots  in  both  cases 
(Figures  k   and  5)  were  used  with  the  eddy  viscosity  model 

vt  -  JT~  i  (78) 

to  test  the  ability  of  this  model  to  reproduce  experimental  results  for 
wakes  in  constant  pressure  gradients.   The  comparison  of  the  analytical 
results  with  experiment  shown  in  Figure  5  verifies  that,  with  some  further 
study  of  the  proper  values  of  the  empirical  constants,  the  approach  out- 
lined here  presents  a  consistent  method  for  testing  wake-type  eddy 
viscosity  models  with  streamwise  pressure  gradients  and  predicting 
untested  physical  situations  with  reliable  eddy  viscosity  models. 

Although  the  analytical  form  of  the  solution  has  been  obtained, 
often  the  initial  conditions  are  discrete  and  do  not  satisfactorily  fit 
any  known  analytic  functions.   In  these  cases  a  numerical  solution  of  the 
governing  equations  can  be  easily  obtained  and  the  computer  program 
needed  to  perform  such  a  computation  is  listed  and  explained  in  the 
Appendix. 

CONCLUSIONS 
The  equations  of  momentum  and  turbulent  kinetic  energy  appropriate 
to  free  turbulent  mixing  have  been  developed.  The  resultant  equations 
have  been  transformed  into  a  plane  which  is  independent  of  the  eddy 
viscosity  model  and  have  been  analytically  solved  by  a  perturbation 
technique.  The  solution  depends  upon  prior  knowledge  of  u  (cp)  and  i(qp) 
and  to  specify  these,  an  approximation  for  cp(x)  must  be  available  in 
order  to  study  a  flow  field  which  is  specified  a  priori.  For  this 
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analysis,  the  wake  model  given  by  Equation  (52)  proved  to  be  satisfactory. 
If,  on  the  other  hand,  only  general  classes  of  flow  fields  are  of 
interest,  u  (cp)  and  4(cp)  may  be  specified  independent  of  a  specific 
physical  problem  with  the  resultant  re -transformation  determining  the 
flow  field  which  they  implied.  The  results  of  the  analysis  have  been 
compared  to  experimental  wake  data  in  constant  adverse  and  favorable 
pressure  gradients.  The  results  of  that  comparison  indicate  that  the 
analysis  supplies  a  satisfactory  method  for  obtaining  analytical  solu- 
tions to  wake  problems  in  freestream  pressure  gradients. 
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FIG.    1 
INITIAL  CONDITIONS   FOR   SAMPLE  CALCULATIONS. 
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SOLUTION  FIELD  FOR  FAVORABLE  PRESSURE 
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C  APPENDIX 

C 

C  THE  FCLLOWING  CCMPUTER  PROGRAM  IS  A  FORTRAN  IV  CODED  FINITE  DIFFERENCE 

C  CALCULATION  OF  THE  ZEROTH  AND  FIRST  ORDER  SOLUTIONS  TO  THE  EQUATIONS 

C  CF  MOMENTUM  AND  TURBULENT  KINETIC  ENERGY  APPLIED  TO  THE  PROBLEM  OF 

C  FREE  TURBULENT  MIXING  IN  FAVORABLE  AND  ADVERSE  PRESSURE  GRADIENTS 

C 

C  TFF  FOLLOWING  PARAMETERS  MUST  BE  INPUT  AT  THE  BEGINNING  OF  THE  DECK 

C  UINF  U(1,J)  EKK1.J)  SIGMA  C  CI  AL  AO  ALL  MY  MX  DELY  DELX  Dl  AKK 

C 

C  *  * 

C  *  UINF  IS  THE  INLET  FREESTREAM  VELOCITY      * 

C  *  * 

C  *  UINF  ALSO  SERVES  AS  THE  PERTURBATION  UO     * 

PA  4e 

C.  *  SIGMA  IS  THE  CONSTANT  IN  ENERGY  DIFFUSION    * 

r  *  if. 

C  *  C  IS  THE  CONSTANT  IN  THE  DISSIPATION  TERM    * 

C  *  * 

C  *  CI  IS  THE  CONSTANT  IN  THE  EDDY  VISCOSITY     * 

r  *  * 

C  *  AL  IS  THE  AREA  OF  THE  TEST  SECTION  AT  X=L    * 

C  *  AO  IS  THE  AREA  OF  THE  TEST  SECTION  AT  X  =  0    * 

C  *  * 

C  *  ALL  IS  THE  LENGTH  OF  THE  TEST  SECTION      * 

C  ~  * 

C  *  MY  IS  THE  NUMBER  OF  LATERAL  POINTS  CALCULATED  * 

C  *  * 

C  *  MX  IS  THE  NUMBER  OF  AXIAL  POINTS  CALCULATED   * 

C  *  * 

C  *  DELY  IS  THE  LATERAL  GRID  SPACING        * 

C  *  * 

C  *  Dl  IS  THE  CONSTANT  IN  THE  STEP  SIZE       * 

C  *  * 

C  *  DELX  IS  THE  AXIAL  SPACING  TO  RETRANSFORM    * 

C  it*  •              ♦ 

C  *  AKK  IS  THE  CONSTANT  IN  THE  APPROXIMATE  PHI(X)  * 

C  *  * 

C  *  MUMMY2  IS  AN  INDEX  FOP  A  PARTICULAR  PSI     * 

C  *  * 


c 


DIMENSION  CHI1 ( 200 , 81 ) , CHI2( 200, 8 1 ) , AK 1 ( 200, 81) , AK2( 200 , 81) 
DIMENSION  THETAC200)  ,DELSTR(200) ,SHAPE (200 ) , PS  I E (200 ) , ELL ( 200 ) 
DIMENSION  PHI (200) ,P2(200) ,P4(200) ,X(200) , FF( 2  00) ,GG (200 ) 
DIMENSION  APS  I  (81) , AU ( 81 ) , AK ( 81 ) , Y( 81 ) , U ( 81 ) , URATI0( 81) ,EK1( 81) 
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DIMENSION    PSK81  ) 

DO    1    1=1, MX 
HO    1    J=1,MY 

CHIK  I,J)=0.0 
CHI2C  I, J)  =  0.0 
AKKI  ♦  J)  =0.0 
AK2U  ,  J)=0.0 

THETAd  )=0.0 
CELSTRt  I  )  =  0.0 
SHAPEC  I  J  =0.0 
PSIEU  )=0.0 
ELL(I )  =  0.0 
PHICI )=0.0 
PSK  I  )=0.0 
P?( I )=0.0 
P4(I )=0.0 
X( I )=0.0 
FP(  I  )  =  0.0 
GGCI  )=0.0 

APSK  J)=0.0 
AU( J)=0.0 
AK( J)=0.0 
Y( J)=0.0 
U( J) =0. 3 
URATIOf J)=0.0 
EK1 ( J)=0.0 

1    CONTINUE 

THE    FCLLOWING    METHOD    OF    INPUTING    U    AND       EK1    USES    A    "STANDARD"    SHAPE 
FOR    ROTH    THE    VELOCITY    AND    TURBULENT    KINETIC    ENERGY    PROFILES    FOR    USE 
IN    TEST     CASES.     IF    THIS     INPUT     IS     DESIRED    ADDITIONAL     INPUTS    ARE    NEEDED. 
UCL-THE    CENTERLINE     VFLQCITY    AT    THE     INLET,     EKCL-THE    TURBULENT     KINETIC 
FMFRGY    CN    THE    CENTERLINE    AT    THE     INLET,    AND    EKMAX-THE    MAXIMUM    TURBULENT 
KINFTIC    ENERGY    AT    THE     INLET.    THESE    VALUES    SCALE    THE    SHAPE    TO    PROVIDE 
SUITABLE     INITIAL    CONDITICNS    FOR    USE     IN    EVALUATING    THE    CONSTANTS     AND 
SPECIFYING    A    STCP    SIZE    WHICH    WILL    BE    BOTH    STABLE    AND    ACCURATE 

U( 1)=100.0 
U(2)=99.0 
U(3)=98.0 
UI4)=96. 0 
U(5)=92.0 
U(6) =88.0 


30 


c 
c 

c 


U(7) 
UI8) 

u(<n 
U(  10 

U(ll 

U(12 
U(  13 
U(14 
U(15 
U(  16 
U(  17 
U(18 
U(  19 
U(20 
U(21 
U(22 
U(23 
U(24 
U(25 
U(26 
U(27 
U(28 


=  81.0 

=  73.0 

=68.0 

)=59.0 

)=52.0 

)=48.0 

)=41.3 

)=35.0 

)=30.0 

)=26.D 

)=21.0 

J=18.0 

)=14.0 

)=11.0 

)=8.0 

)=5.0 

)=3.0 

)=2.0 

)  =  1.0 

)=0.5 

)=0.25 

)  =  0.0 


ALPHA=1.3~(UCL/UINE) 


DD    2    J=l,51 
U( J) =UINF*(1.0 
2    CONTINUE 


^LPhA"U  (J)/100.0  ) 


FKK 
EKK 
EKK 
EK1( 
EKK 
EKK 
EKK 
EKK 
EKK 
EKK 
EKK 
EKK 
EK1( 
EKK 
EKK 
EKK 
EKK 
EKK 
EKK 


1) =80.3 

2)  =  81.0 

3)=82.0 

4)  =84.0 

5)=87.0 

6)=90.0 

7)=95.0 

8)=98.0 

9J=100.0 

10)=98.0 

11)=92.0 

12)=83.0 

131=75.0 

14)=67.0 

15)=60.0 

16)=53.0 

17)=48.0 

18)=41.0 

19) =36.0 
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c 


EK1(20)=32.0 

EKK21  )=26.0 

FK1(22)=21.0 

EK1C 23)=17.0 

EK1(24)=13.0 

EKK25)  =10.0 

EK1(26)=8.0 

EKK27)=6.0 

EK1<28)=4.0 

EKK  29)  =  2.0 

EK1(30)=1.0 

EK1(31)=0.25 

EK1(32)=0.0 

GAMMA=5.0*(1.0-(EKCL/EKMAX) ) 

DO    344    J=l,9 

EKK  J)  =  EKMAXM1.0-(1.0-EK1  (  J  )/  100  .0  )*GAMMA  ) 

344  CONTINUE 

DO    345    J=10,51 

EK1( J)=CKMAX*EK1 (J) /l 0  3.0 

345  CONTINUE 
C 

UINF2=UINF**2 

C=C/C1 

BETA={  (AO/AL )**2-l .0  )/ ALL 
C 
C  SET    UP     INITIAL    Y    FIELD 

Yd  )=0.0 

DO    50    J=2,MY 

Y(J)=FLOAT( J-l )*DELY 

50  CONTINUE 
C 

C  SET    UP     INITIAL     PS  I     FIELD 

PSI (1 )=0.0 

PSI (2)  =  (U(  1)+U(2))*0.5*DELY 
DO    51    K=3,MY 
K1=K-1 

A=0.5*(U( 1)+U(K) ) 
B  =  0.0 

DC    52    L=2,K1 
B=U(L)+B 
52     CONTINUE 

PSI  (K)=DELY*  (A  +  B  ) 

51  CONTINUF 


NY1=MY-1 
DELPSI=PSI(MY) /MY1 
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C      SET  UP  UNIFORM  PSI  FIELD 
APSI (1) =0.0 
DO  53  J=2,MY 
APSK  J)=FLOAT{  J-1)*DFLPSI 

53  CGNTINUE 
f. 

WRITE(6,995) 

995  F0RMAT(1H1,40X, *****  TABULATED  INPUT  VARIABLES  IN  PHYSICAL  CCORDIN 
1ATES  ****• ) 

C 

WRI"<"c  (6,770) 
770  FGRMAK 1H0.25X,1 VELOCITY' ,10X, 'TURBULENT  KINETIC  ENERGY ', 10X ,' L AT E 
1RAL  POSITION' ,10X, 'STREAM  FUNCTION') 
C 

00  880  J=1,MY 

WRITE (6,996)  U(J),EK1(J),Y(J),APSI(J) 

996  FCRMAT(1H0,2  5X,E15.7,10X,E15.7,10X,E15.7,10X,E15.7) 
880  CONTINUF 

C 

C  REDISTRIBUTE    THE    INPUT    VALUES     IN    A    UNIFORM    PSI    FIELD 

AU(1)=U(1) 

AK(1 )=EK1(1 ) 

DO    54    J=2,VY 

PCS=APSI(J) 

AU  (  J  )  =  P  I  F2  (  POS  ,  P  S  U MY ,  U  ) 

AK( J)=PIF2(P0S,PSI , MY,EK1) 

54  CONTINUE 
C 

C  SET    THE    VALUES    FOR    THE    DEPENDENT    VARIABLES    OF    THE    CALCULATION 

DO    55    J=1,MY 

CHIH1,  J)=UINF  **2-AU(J)*J-2 
AKK  1,  J)=AK(  J) 
PSI(J)=APSI(J) 
U( J)=AU( J) 

55  CONTINUE 
C 

DO    40    J=1,MY 

MUMMY=J 

MUMMY1=J-1 

IFCCHI1 (1, J) .LT.0.1 )    GO    TO   41 

40  CONTINUE 

41  CONTINUE 
C 

PSIEll)=PSI( MUMMY) 

UF    =U(MUMMY) 

UE2=UE=**2 
C 
C      SET  UP  NEW  Y  FIELD  BASED  ON  UNIFROM  PSI  FIELD 
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c 
r. 
c 


Y(1)=0.0 

Y(2)=((U(1)+U(2) )/(2.0*U(l )*U(2)  ) KDELPSI 

DG    5fc    K=3,MY 

K1=K-1 

A  =  0.5*(  (1.0/UC1)  )+(1.0/U(K)) ) 
.  B=0.0 

DO    57    L=2,K1 

B=(1.0/U(L) )+B 
57    CONTINUE 

Y(K)=DELPSI*(A+B) 
56    CONTINUE 

ELLC1 )=Y(MUMMY  ) 


DO    43    J=1,MY 

URATIO( J)=SORT(1.0-CHIl( 1, JJ/UE2) 

43  CONTINUE 
C 

C      CALCULATE  DISPLACEMENT  THICKNESS 

A=0.5*{ ( 1.0/URATIO(1))+ (1.0/ URATIO (MUMMY) )-2.0) 

B  =  0.0 

CO    44    J=2,MUMMY1 

B=( ( 1.3/URATI0( J))-1.0)+B 

44  CONTINUE 
DFLSTR(l)=DELPSI*(A+8)/UE 

C 

C  CALCULATE    MOMENTUM    THICKNESS 

A=0.5*(2.0-URATIC(1 )-URAT 10 ( MUMMY ) ) 

R  =  3.0 

DO    45    J=2,MUMMY1 

B=(1.0-UPATI0< J) )+B 

45  CONTINUE 
THETAU  )=DELPSI*(A+B)/UE 


THETOTHETM 1) 

SHAPF(1)=DELSTRC1 )/THETA(l ) 

WRITE(6,997) 
997    FORMAT ( 1H1 f40Xf *****    TABULATED     INPUT    VARIABLES    IN    STREAM    FUNCTION 
1CC10RDINATES    **"**'  ) 

WRITE (6,666)     PHI ( 1 )  ,TH ET A ( 1 ) , S HAP E( 1 ), DELS TR ( 1 ) , ELL ( 1) , PS  IE ( 1 ) 

666    FORMAT( 1H0,2X, 'PHI     =• , E 10. 4 , 2X , ' THETA    = ' , E 10. 4 ,2X , ' SHAP E    FACTOR    =  • 
1 ,P10.4,2X,'DELSTAR    = • , E 10.4, 2X, • MACROSCAL E    = ■ . E 10. 4, 2X, » MA SS    FLUX 
2=« ,E10.4) 

WP.ITEC6.771) 
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771  FORMAT! 1H0,20X, »CHI',15X, 'VELOCITY' ,10X, 'TURBULENT  KINETIC  ENERGY' 
ltlOXt 'LATERAL  POSI T ICN' , 10X ,' STREAM  FUNCTICN') 
C 

CO  881  J=1,MY 

WRITE  (6, 99  8)  CHIHl  ,J)  ,  AU ( J)  , AK1 < 1 , J ) , Y ( J ) , PS  I  ( J ) 
998  FORM AT ( 1 HO, 10X, El 5. 7,1  OX, El 5. 7,1 OX, El  5. 7,1  OX,  El  5. 7,1  OX, El  5. 7) 
881  CONTINUE 
C 
C      SET  STABLE  STEP  SIZE  FOR  DELPHI 

DELPHI=D1*(DELPSI**2) 
C 

PHI(1)=0.0 
CO  5  1=2, MX 

PHK  I)=FLOAT(I-l)*DELPHI 
5    CONTINUE 


AMEGA=0.5*BETA/(AKK*THETO* (UINF**2) ) 
C 

CO    7    1=1, MX 

P2(I)  =  1.0/(  <1.3+AMEGA*PHIU)  )**2) 

P4(I)=1.0/( { 1.0+AMEGA*PHI(I) )**4) 
7    CONTINUE 
C 

MX1=MX-1 
C 
C  START    MAIN    CALCULATICN    LOOP 

DO    20    1=1, MX1 
C 

F=-C/((PSIE( I)+UINF*THET0*P2 (I ) )**2) 
C 

DO    10    J=2,MY1 

DERIV=(CHI1(I,  J  +  D+CHIKI  ,  J-l )  -2.0*CHI  1  (  I ,  J  )  )  /  (  DEL  PS  1**2  ) 

CHIKI  +  1,  J)=CHI1(  I,  J)+DELPHI*DERIV 
10    CONTINUE 


C 

c 


DERIV=2.0*(CHI 1( I,2)-CHI1( I, 1) ) / ( DELPSI**2 ) 
CHI1(I+1,1)=CHI1(I,1)+DELPHI*DERIV 

DO    11    J=2,MY1 

G=-0.5*P2( I)*CHI1(I,J)*(CHI1(I,J+1)-2.0*CHI1(  I , J )+CHIl(  I  ,  J-l )  )  / 
1(DELPSI**2) 
DERIV=(CHI2( I,J+1)-2.0*CHI2( I , J ) +CHI 2( I , J-l ) ) /( DE LP SI** 2 ) 
CHI2(I+1,J)=CHI2(I,J)+DELPHI*(DERIV+G) 
11    CONTINUE 

G=-0.5*P2(I  )*CHU(I  ,1)*2.0*(CHI1(I,2)-CHI1(  1 ,  1 )  )/  (  DELPS  1**2  i 
DERIV=2.0*(CHI 2( I,2)-CHI2( I ,1) )/(DELPSI**2) 
CH 12 (1+1,1 )=CHI2( I,11+DELPHI*{DERIV+G) 
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c 
c 


DO    12    J=2,MY1 

G=0.25*P2(  I)*(  (  (CHIKI  ,J+1)-CHI1(I  ,  J-l ) ) / ( 2.0*DELPS I*U I NF)  )**2 ) 

CERIV=(AK1(  I,J+1)-2.0*AK1(  I , J ) +AK 1 ( I , J- 1 ) ) / ( DELP S 1**2 ) 

AKKI  +  1,  J)  =AK1  (  I,  J)+DELPHI*{  ( DER I V/S I GMA ) + F* AK1 {  I,  J  )+G) 

12  CONTINUE 

DERIV=2.0*(AK1 ( I ,2 )-AKl (1,1)  )/ ( DEL PS  1**2 ) 

AKKI  +  1,  1)=AK1(  I,1)+DELPHI*(  (  DER  I  V/S  I  GMA  )  +F*  AK1  (  I  ,  1  )  ) 

DO    13    J=2,MY1 

G1A  =  CHIKI,J)*(  AKK  I  ,  J  +  l  )-2.  0*AK1  (  I  ,  JJ  +  AKKI  ,J-1)  )  /  (  DEL  PS  1**2  ) 

G1B=(CHU(  I,J  +  1)-CHIK  I,J-1)  )*(AKl(I  ,J+1)-AK1(I,J-1)  )/(4.0* 
1(DELPSI**2)) 

G1=-0.5*P2(I )*(G1A+G1B ) 

G2=0.12  5*P4(  I)*CHI1(I,  J)*(  (  (  CHIKI, J+l) -CHIK  I , J- 1  )  ) /( 2 . 0*DE LPSI * 
1UINF)  )**2) 

G3=0.5*P2(  I)*  (CHIK  I,J  +  1)-CHIK  I  ,  J-l )  )* (  CH  I  2  (  I  ,J+1)-CHI2(I  ,J-1)  )/ 
1 (4.0*(DELPSI**2)*(UINF**2) ) 

G4=-0.5*C*P4(I  )*CHIKI,J)*AK1(I,J)/((PSIE(  I  )+UI  NF*THET0*P2(  I ) )**2) 

G=G1+G2+G3+G4 

DERIV=(AK2(I,J+1)-2.3*AK2( I, J)+AK2( I. J-l) ) / ( DELPS 1**2 ) 

AK2( 1+1, J)=AK2( I, J )+DELPHI*( ( DERI V/S I GMA ) +F* AK2 ( I ,J)+G) 

13  CONTINUE 

G1A=CHI  1(1,  1)*2.0*(AK1(I  ,2)-AKKI  ,1)  )  /  (  DEL  PS  1**2  ) 

G1B=0.0 

G1=-3.5*P2( I)*(G1A+G1B) 

G2=0.0 

G3=0.0 

G4=-0.5*C*P4(I )*CHI  1(1  ,1)*AKK 1,1 )/( (PSIE(  I  )+UINF*THETO*P2( I  ))**2) 

G=G1+G2+G3+G4 

CERIV  =  2.0*(AK2( I,2)-AK2( 1 , 1 )  ) / ( DEL  PS  1**2 ) 

AK2(  I  +  1,1)=AK2(I  ,  D+DELPHI*  (  (  DERI  V/S  I  GMA)  +  F*AK2  (  1,1  )+G) 
C  END    OF    MAIN    CALCULATION    LOOP 

C 

DO    75    J=1,MY 

U( J)  =  SORT(UINF2/P2(  I  )-CHI K 1  +  1 , J ) ) 

AU( J)=SORT(UINF2/P2 ( I )-( CHI 1 (  1  +  1 , J ) +CHI2 ( I +1, J )/U INF2) ) 
75    CONTINUE 
C 

DO    61    J=1,MY 

MUMMY=J 

MUMMY1-J-1 

IF(CHI1(I+1,J) .LT.0.1 )    GO   TO    70 
61    CONTINUE 
C 

70    CONTINUE 
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c 
c 


MUMMY2=MUMMY/2 

PSIE  (1+1 )  =  PSI (MUMMY) 


.    UC=U(MUMMY) 

UE2=UE**2 
C 
C  CALCULATE    NEW    Y    FIELD 

Y(1)=0.0 

Y(2)=0.5*DELPSI*(1.0/U(1)+1.0/U(2) ) 

DO    59    K=3,MY 

K1=K-1 

A=0.5*(  1.0/U(1)+1.0/U(K)  ) 

B  =  0.0 

DO    60    L=2,K1 

B=1.0/U(L)+B 
60    CONTINUE 

Y(K)=DELPSI* (A+B) 
59    CONTINUE 


ELL(I+1)=Y(MUMMY) 

DO    62    J=1,MY 

URATIOC J)=SQRT(1.D-CHI1(I+1, J )/UE2) 

APSI(J)=SQRT(1.0-(CHI1(I+1,J)+CHI2(I+1,J)/UINF2)/UE2) 

EKK  JJ  =  AK1(  1+1,  JJ+AK2(  1+1, J)/UINF2 

62  CONTINUE 

A=0.5*( (1.0/URATI0{1)  )  +  ( 1.0/URAT 10 ( MUMMY) )-2.  0) 

B  =  0.0 

DO  63  J=2,MUMMY1 

B  =  (  (1.0/URAT IOC  J) )-1.0)+B 

63  CONTINUE 

DELSTR( I+1)=DELPSI*(A+B)/UE 

A=0.5*(2.0-URATI0( 1 J-URATIO ( MUMMY) ) 

B  =  0.0 

DO    64    J=2,MUMMY1 

B=( 1.0-URATIO( J) )+B 

64  CONTINUE 
THFTA(I+1)=DELPSI*(A+B)/UE 

SHAPE (  I  +  1)=DELSTR( 1  +  1) /THETAd+1) 

WRITE(6,991) 
991    FCRMAT(lHlf40Xf •****    INTEGRATED    PROFILE    PARAMETERS    ****•) 


37 


WRITE (6,992)  PHI (1+1)  ,THETA( I+1),SHAPE(  I  +  1 ) , DELSTR ( 1  +  1 ) , ELL ( I +  1 ) , 
1PSIEU  +  1) 

992  FQPMAT(1H0,2X, 'PHI  = • , E 10 . 4, 2X ,  'THET A  = • , E 10. 4, 2X, ■ SHAP E  FACTOR  =• 
1,E10.4,2X,»DELSTAR  =• , E10.4 , 2X , ' MACROSCAL E  = ' , E10.4, 2X, 'MASS  FLUX 
2=»,E10.4) 

WRITE(6,999) 
999  FCRMAT( 1H0,40X, •****  CALCULATED  OUTPUT  VARIABLES  ****■ i 

WRITE(6,772) 
772  FORMAT ( lHOt  »J« ,8X,  'U/UE • ,8X,'UT/UE'  , 8X,'U«  » 8X, " UT' , 8X t* Y"  ,8X,»PSI» 
1  ,8X,  'CHIl'  ,8X,  •CHI2S8X,  'AK1S8X,  •  AK 2  '  ,  8X,  »AKT«  ) 

DO  883  J=1,MY 

WRIT  EC 6, 990)  J , URAT  10 ( J ) , APS  I  ( J  ) , U ( J ),AU(J),Y(J),PSI(J)tCHI1(I+1,J 
D,CHI2(I  +  ltJ)»AKl(I+lf  J)  ,AK2(I+1, J)  ,  EKK  J) 
9  90  FORMAT (1 HO, I2,2X,F5.3,2X,F5.3,2X,F7.3t2X,F7.3t2X,F5.3,2X,F7.3,2X,E 

112.5,2X,E12.5,2X,E12.5,2X,E12.5,2X,E12.5) 
883  CONTINUE 
ICALC=I 
IF(U(1) .LT.5.0)  GO  TO  21 

20  CONTINUE 

21  CONTINUE 
WRITE(6,993) 

993  FORMAT( 1H1,1****  SUMMARY  OF  INTEGRATED  QUANTITIES  ****•) 

DO  874  I=1,ICALC 

WRITE (6,994)  THE TA ( I )  , DELSTR ( I )  , SHAPE ( I )  ,  PS  I  E  (  I ) , ELL ( I ) , I 
C94  FORMATClHOt 5E20.8, 13) 

874  CCNTINUF 

X( 1)=0.0 
DC    869     I =2, MX 
X( I)=FL0AT(I-1)*DELX 
8  69     CONTINUE 

FF(1)=0.0 
DO    870    1=2, MX 

FF(I)=2.0*C1*UINF*( (1.0+BETA*XC I ) )**1 .5-1 .0 )/ ( 3.0*BETA) 
870    CONTINUE 

GG(1 )  =  3.0 

GG(2)=0.5*(     1. 0/(SORT(AKK 1,MUMMY2) )*ELL(1) )+    1.0/ (SORT (AK1(2,MUMM 
1Y2))*ELL(2)) )*DELPHI 
DO    871     I=3,ICALC 
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11=1-1 

A  =  0.5*(  1.0/ (  SORT  (AKK  1  ,MUMMY2  ) )  *E  LLC  1)  )+ 1.0/ (SORT  ( AK1 ( I , MUMMY2 )  )* 
1ELL(I)  )  ) 

3  =  0.0 

DO  872  K=1,I1 

B=1.0/(  SORT  (AKK  I,MUMMY2)  )*ELL(K)  )+B 
872    CONTINUE 

GG(I  )  =  DELPHI*(A  +  B) 
871     CONTINUE 

DO    877    I=ICALC,MX 

GG(  I  )=0.0 
877    CONTINUE 


850 


848 


873 
849 


WRITE(6,850)    CI 

FORMAT (lHl,30Xt«RETRANSF0RMATI0N    RESULTS    FOR    CI    =',E10.3) 

WRITE(6,848) 

FORMAT (1 HO, 30X  , «F« , 20X,«G« ,  20X, 'PHI' ,  20X, »X' ) 

DO    849    1=1, MX 

WRITE (6, 873)    FF(I),GG(I),PHI(I),X(I) 

FORMAT! 1H0,10X,4E2  0.8) 

CONTINUE 


END 


FUNCTION    PIF2     (X ,XL 1ST , N, FL 1ST ) 
FUNCTION    PIF2    IS    A    SECOND    ORDER    LOOKUP    FUNCTION 
DIMENSION    XLIST     (100),     FLIST    (100) 
BLIF     (P,0,R,S,T)     =     UQ-P)*(S-T)/(R-Q)+S) 
IF    (X-XLIST(N) )     2,1,1 
1    I    =   NtI 
GO    TO    5 

IF(X-XLIST( 1))     4,4,6 
I    =    1 
K    =    1 
GO    TO    30 
K    =   2 

DO  8  I  =  1,N 
IF  (X-XLIST( I) )  9,9,8 
CONTINUE 
I=N 

I  =  1-1 

BLIF1  =  BLIF(X,XLIST(I ) ,XL 1ST ( 1+1 ) , FL I  ST ( I  ),FLIST(I  +  1) ) 
IF  (K-l)  11,11,12 
PIF2  =  BLIF1 
RETURN 


2 
4 
5 

6 

7 


9 
30 
10 
11 
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12 
13 
14 
15 

16 
17 

18 


IF((I+2)-N)  13,13,16 

IF  (CI-l)-l)  15,14,14 

IF(ABS(XLIST(I-l)-X)-ABS(XLIST(I+2)-X)) 16,15,15 

L  =  1+2 

GO  TO  17 

L  =  1-1 

BLIF2  =  BLIF  (  X  ,  XL  I  ST  (  I  )  ,XL  I  ST  (  L)  ,  FL  1ST  (  I  )  ,  Fl_  I  ST  (  L)  ) 

PIF2  =  BLIF  (X,XLI STCI+1) ,XLIST(L) ,BLIF1 ,BLIF2) 

RETURN 

END 


ko 
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